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Abstract: We study the evolution of density perturbations for a class of f(R) models 
which closely mimic ACDM background cosmology. Using the quasi-static approximation, 
and the fact that these models are equivalent to scalar-tensor gravity, we write the modified 
Friedmann and cosmological perturbation equations in terms of the mass M of the scalar 
field. Using the perturbation equations, we then derive an analytic expression for the 
growth parameter 7 in terms of M, and use our result to reconstruct the linear matter 
power spectrum. We find that the power spectrum at z ~ is characterized by a tilt 
relative to its General Relativistic form, with increased power on small scales. We discuss 
how one has to modify the standard, constant 7 prescription in order to study structure 
formation for this class of models. Since 7 is now scale and time dependent, both the 
amplitude and transfer function associated with the linear matter power spectrum will be 
modified. We suggest a simple parameterization for the mass of the scalar field, which 
allows us to calculate the matter power spectrum for a broad class of f(R) models. 
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1. Introduction 

Cosmological data, arising from anisotropics of the cosmic microwave background [1, 2], 
large scale structure [3] and measurements of type la supernovae [4, 5], indicate that the 
expansion of the Universe is currently accelerating. The simplest approach to modeling this 
epoch is to postulate the existence of a very small but non-zero vacuum energy, which can 
dominate at late times and accelerate the expansion by virtue of having negative pressure 
(specifically, an equation of state parameter w = —1.) However, the extreme fine tuning 
implicit in such a model has driven a search for alternative dark energy candidates, where 
the acceleration is attributed to one or more dynamical fields [6], (alternative approaches 
include higher dimensional physics [8, 9, 10] or the backreaction of inhomogeneities on an 
isotropic, averaged metric; see [11] for a recent discussion.) 

Rather than introducing additional matter content to the Universe, an alternative 
approach to modeling the current epoch is to postulate that at cosmological distance scales, 
gravity deviates from its standard General Relativistic description. This can be achieved, 
for example, by introducing fields which mediate gravity in addition to the standard spin-2 
graviton. In this work we consider the f(R) subclass of scalar-tensor theories, described 
by the action 
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where C m is the Lagrange density of standard matter and radiation, and f(R) is an un- 
specified function of the Ricci scalar 1 . These models contain an additional scalar field 
in the gravitational particle sector, which we will call the scalaron. The scalaron has 
the unique property that its mass depends on the background curvature of spacetime 
M = M(i?background); this fact proves important in evading solar system tests of gravity. 
A large body of literature has been devoted to models of the form (1.1), and we direct the 
reader to [7, 12, 13, 14] and references therein for a detailed recent review. 

In principle one could use an arbitrary f(R) function in the action (1.1), however in 
order to be both observationally and theoretically viable these models must respect a large 
list of consistency requirements. To begin, the f(R) function must satisfy 1 + /r > 
and /rr > (throughout the paper, R subscripts denote derivatives with respect to R) 
in all dynamically accessible regions to cosmology, to evade ghost and other instabilities 
[15, 16, 17, 18, 19]. 

f(R) models must also have a viable Newtonian limit. The conditions required for an 
f(R) function to reproduce Newtonian gravity at a particular curvature scale i?N are 



\f{R N )\^R N |/r|<1 RnIrr^L (1.2) 

These conditions must hold, for example, in the Solar System, where gravity is accurately 
described by weak field Newtonian physics. 

Finally, prospective f(R) models must be able to reproduce a viable cosmology [20, 
21, 22], with a period of radiation domination followed by a standard matter era where 
the scale factor evolves approximately as a(t) ~ i 2//3 . It has been shown that only a very 
restricted subclass of f(R) models can reproduce an expansion history consistent with 
observations, as the scale factor typically evolves as a(t) oc i 1 / 2 during matter domination 
[20, 21]. 

In addition to the above theoretical considerations, Solar system and laboratory tests 
of gravity provide stringent observational constraints on f(R) models [24], [23], [25]. For 
example the Cassini probe measurement of the parameterized post Newtonian (PPN) pa- 
rameter 7ppn in the solar system [26] can be translated into a constraint on |/jj(-R g )| [24] 

\f R (R g )\ < 4.9 x 10- 11 , (1.3) 

where R g ~ 87rGp g and p g ~ 10~ 24 g cm -3 is the typical galactic background energy 
density. The condition (1.3) might appear to be a particularly strong condition on f(R) 
models (/# = is the General Relativistic limit), however the 'background' dependence 
of the mass of the scalar field allows us to construct f(R) models that easily evade this 
bound. Viable models are constructed such that the (squared) mass of the scalar field is 

1 The sign conventions in this paper are: the metric signature (-+++), the curvature tensor R" iipv = 
dpT'liv ~ •■•! Rjxv = FFuavi so that the Ricci scalar R = Rfc > for the de Sitter space-time and the 
matter-dominated cosmological epoch. In addition, we set c = 1 throughout. 
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much larger than the local curvature, M 2 (R g ) 3> R g . In this case, the scalar field will not 
propagate on macroscopic distance scales, and deviations from GR will be suppressed in 
the solar system. 

More recently it has been argued that stringent observational constraints on f(R) 
gravity arise by considering structure formation, as the evolution of density perturbations 
will be significantly modified by the scalar field. The most commonly studied model in 
this regard is the 'minimal' f(R) model [27], so called because it has an expansion history 
that is identically ACDM by construction (although see [30, 31]), and possesses only one 
additional parameter as compared to the standard cosmological model. This additional 
parameter is taken to be Bq, defined as Bo = B(a = 1) where the function B(a) is given 
by 



£ = T^r^r> (1-4) 



/rr R'H 
1 + /r H' 

where primes denote derivatives with respect to ln(a). In a series of papers [28] [29] [32], 
it has been concluded that current cosmological data sets impose [32] Bq < 1.1 x 1CT 3 at 
95% C.L 2 (B = is the General Relativistic limit.) 

Observational tests of f(R) models place lower limits on the mass of the scalar field. 
If the mass is much larger than the background curvature, the scalaron will essentially be 
non-dynamical, and this class of models will mimic General Relativity. In the solar system, 
the constraint (1.3) implies that the scalaron mass must be very large; we would like to 
consider the extent to which the mass is allowed to relax on cosmological scales. It is the 
aim of this work to construct a parameterization of f(R) models, which describes both 
the expansion and growth histories, in terms of the scalar field mass M(a). We derive all 
results in a model independent manner (all quantities will be written in terms of M(a)), 
however we will sometimes resort to specific models when we wish to numerically evolve 
the field equations. We will use the following two functions 



/(fl) = -^+2 9 «-W/-'> + ^, (1.5) 

/(*) = -^ + A^(^) 2 " + ^, (1.0) 

where e = R vac /2g(b + log [2 cosh [6]]), g, n, A, b are model parameters, we will take Mq to be 
an inflationary mass scale and R va , c = 12HqCIa is the curvature associated with the vacuum. 
These two functional forms, which will be referred to as the exponential and power law 
models respectively, represent expansions around General Relativity that are commonly 
studied in the literature [24], [38], [37], [39], [40] [41], [42]. We stress that these functions are 
not globally viable f(R) models, however they are sufficient for our purposes in the sense 
that they are free from instabilities, and reproduce the standard cosmology, over redshifts 
z = (10 3 ,0). 

2 see also [33], [34], [35], [36] for constraints on f(R) models 
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The paper will proceed as follows. We begin by studying the growth parameter 7. 
In section 2.1 we discuss 7 as a parameter in General Relativity. In section 2.2 and 2.3 
we consider the field equations for f(R) models, and how they can be simplified using 
the 'quasi-static' approximation. We use this approximation scheme to construct a model 
independent expression for 7 in terms of the mass of the scalar field M(a). In section 3 
we introduce a simple yet representative parameterization of M(a) and discuss how the 
matter power spectrum depends on this function. We end with a discussion on the use of 
7 as a parameterization for growth, and how our approach relates to existing work in the 
literature. 



2. The growth parameter 



We begin by studying the growth parameter 7 [43, 44, 45, 46, 47], defined as 



7 = 



log[dlog[(5 m ]/dlog[a]] 
log[O m ] 



(2.1) 



where l rrf subscripts denote the total matter component, 5 m = 8p m / p m and J7 m is defined 
in the usual manner, 



fir, 



8vrG/9 ri 



(2.2) 



2.1 The growth parameter in General Relativity 



We first briefly re-derive the value of 7 in a standard ACDM cosmology; 7 ~ 6/11 over the 
redshift range z ~ (0, 10). To obtain this result, consider the evolution of density perturba- 
tions in a cosmological background containing two perfect fluids; pressureless matter (which 
will dominate for z ~ (1, 10)) and a subdominant component with unspecified equation of 
state P = wp (throughout this paper we will make the simplifying assumption that the 
Universe is spatially flat, J\o = 0). Considering only the evolution of matter perturbations 
in this epoch, <5 m evolves according to 



a 2 C + 



3-^ m -^(l + ^)(l-O m ) 



aSL 



3 i 

2 



0, 



(2.3) 



where throughout the paper primes denote differentiation with respect to a. By differenti- 
ating the definition of 7 (2.1) with respect to a, and using (2.3) to remove second derivatives 
of S m , we obtain an equation describing the evolution of 7, 



I '-'--in / •■'--)!) ,, \ --11.) --ju 3. s^m 

7 = 07 i„„ro i I o— + ( - 2 ) — — + 2^ 1 + ~ 




3w 7 
a log [SI 

To solve (2.4), we use the fact that for z = (1, 10), we can expand f2 m as 



(2.4) 
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n m = 1 - ^a~^ + O ( J|pa- 6 ^ , (2.5) 

where fi m o and f2#o are, respectively, the fractional densities of the matter and subdominant 
fluid at the present time. Expanding equation (2.4), we obtain a series solution 

3(1 -w) „ ( ^wO -3tD 



' 5-6^ + °{^ a I' (2 ' 6 ' 

At first order in the expansion 7 is sensitive only to the equation of state w of the sub- 
dominant energy component. For z = (1,10), we take w = — 1 and find 7 ~ 6/11 as 
expected. 

At low redshift z < 1 , we can no longer use the expansion a <C 1 to derive an expression 
for 7. However, it is straightforward to show that for a ACDM cosmology, 7 will asymptote 
to a value 7—7-2/3 for a > 1. This asymptotic behaviour is logarithmic with respect to 
the scale factor, and hence we can assume that 7 will not deviate significantly from the 
value 7 ~ 6/11 for < z < 1. At high redshift z > 10, constructing an analytic form for 7 
is complicated by the fact that we must take into account the radiation component in both 
the perturbation equations and O m (o) (since p r > p\ for z > 10, where p r is the radiation 
energy density.) 

For dark energy models with a constant equation of state, 7 is an excellent param- 
eterization for the growth [46], [4 7]. For redshifts of interest to structure formation, 7 is 
approximately constant, 07' <C 7, and is only weakly dependent on cosmological parame- 
ters such as ftmo. In addition, knowledge of 7 is sufficient to reconstruct the time evolution 
of the matter power spectrum for z < 10. This is achieved by integrating the expression 



S m = exp 



/ 



fH, d In a 



(2.7) 



2.2 The growth parameter in f(R) models 



The above results are well known in General Relativity [43, 44, 45, 46]. We now perform a 
similar calculation for f(R) models. We note that the growth parameter in the context of 
modified gravity has been considered in a number of recent works [49, 50, 51, 52] (see also 
[47] for an earlier treatment of 7 in modified gravity models.) 
At zeroth order in the perturbations, we have 



3H 2 = 87rG( Pm + Pr ) + JR - J - 3Hf R , (2.8) 
-2H = SvrGK + + / R - Hf R , (2.9) 

where dots denote derivatives with respect to time and p m ,Pr are the energy densities 
of matter and radiation respectively. The f(R) terms in (2.8,2.9) are often written as 
components of a perfect fluid, 
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3i7 2 = 87rG(p m + p r + / 9 f ), (2.10) 
-2H = 8wGp m + + 87rG(p f + P f ), (2.11) 

where pt and Pf are effective density and pressure terms 

SnGpi = IkR ~ f - 3Hf R , (2.12) 

SnGPf = - {RfR ~ f) + f K + 2-ff/R. (2.13) 
The equation of state parameter u> f for this fluid is given by 

2f R + 4Hf R -(f R R-f) 

Wi = ■ j (2.14) 

-6F/ R + (f R R -f) 

from which it is clear that the modified gravity terms will mimic dark energy with equation 
of state Wf ~ —1 whenever f R is approximately static; Hf R , f R <C Rf R — f. Although we 
have rewritten the modified gravity terms as pf and Pf, we stress that these quantities are 
functions of H, H, H, and the modified gravity equations are fourth order in derivatives of 
the scale factor. 

We are interested in the behaviour of the density perturbations 5 m in this class of 
models. The equations describing the evolution of scalar perturbations in f(R) models 
have been derived in [54]. In the Newtonian gauge, using sign conventions such that the 
metric potentials are defined as 

ds 2 = -(1 + 24>)dt 2 + a 2 (t)(l - 24>)i ij dx i dx j , (2.15) 

we find [54] 

S m + 2HS m + fcV - 30 - 6F0 = 0, (2.16) 

(1 + MM -</>) + fate = 0, (2.17) 



(i + /i 



2^0 + 6F(0 + #V) 



a 2 



,'k 2 



+ 3/rr(P + H z ) - ^/rr + 3Hf RR e + 



a 2 ' 



-3Hf RR e + /r(6^V + 30) = -8irGp m 5 m , (2.18) 

k 2 ■■ ■ k 2 

e = -12(H + 2H 2 U - 6Hip + 2-^ - 6(0 + 4F0) - 4^,0, (2.19) 

a z a z 

where we have neglected perturbations in the radiation component (this is a valid approx- 
imation at late times.) 

The background and perturbation equations (2.8,2.9,2.16 — 2.19) comprise a com- 
plicated system of coupled, fourth order differential equations, and solving them (even 
numerically) is a highly non-trivial task. However, for the class of so-called 'viable' models 
which can be written as expansions around General Relativity (see for example (1.5, 1.6)), 
we can use the quasi-static approximation [57], [58] to simplify the dynamics. 
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2.3 Field equations in the quasi-static approximation 

In a cosmological context, the quasi-static approximation posits that the mass of the scalar 
field remains larger than the Hubble parameter at all redshifts; H(a)/M(a) <C 1. In this 
case, the scalar field will only be weakly dynamical throughout the cosmological history 
and the background evolution will closely mimic ACDM, with corrections of order H 2 /M 2 . 

To see what happens to the modified Friedmann equation (2.10) in the quasi-static 
limit, we solve the trace of the gravitational field equations 

3/rrR + 3/ RRR (i?) 2 + 9H/ RR R + R - R/ R + 2/ = -SttGT. (2.20) 

For models that reduce to expansions around General Relativity (such as (1.5,1.6)), we 
can expand the function f(R) around a cosmological constant, f(R) = — -R vac /2 + f(R) , 
where / is an unspecified function of the Ricci scalar. Doing so, it has been found [38] that 
equation (2.20) has approximate solution R ~ Rqr + $Rmd + 5R OS c, where -Rgr and SR{ n d 
are given by 



12i^=4^ + ° A ' (2 - 21) 
SRmd = — 3/rr-Rgr — 3/ RRR (.R GR ) 2 — 9H / RR -Rgr + -Rgr/r — 2/. (2.22) 

In (2.22) and unless stated otherwise throughout the paper, all /, /r, /rr, /rrr terms are 
functions of R GR ; f = /(Rgr)- 

5R osc is a component of the Ricci scalar which oscillates with high frequency oj ~ M(a). 
The oscillating component has been discussed in previous works, see for example [38, 55, 56], 
and it has been argued that the energy density associated with these oscillations will decay 
in the early Universe, making SR osc negligible at late times z < 10 3 . We therefore neglect 
<5 i?osc in what follows. In addition, 5i?i nc i satisfies 5R m d <C -Rgr in the quasi-static limit, 
and hence we can also neglect this term. We are left with R ~ -Rgr> which can now be 
substituted into the Friedmann equation. We find 

H 2 n m0 n i0 R G r j / H^ R , - 

Hl-^ + ^ + nA0+ QHl fK ~ 6H 2 ~ G 1^^R/RR, (2-23) 

where i?GR is given by 

^ = ^ + ^? + fW (2.24) 
Hq a 6 a 4 

The equation of state of the /(R) effective fluid can be written as 
w ^ 1 , aHGR(aHGR/ R )' aH^R'^f^ 



3 It is important to stress that there is no true cosmological constant in these models; they generically 
admit Minkowski space as a solution to the vacuum field equations. However, for the region of interest to 
cosmology we can expand the f(R) function around a pseudo constant 7? vac /2. 
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The f(R) terms on the right hand sides of (2.23,2.25) are of order Hq R /M 2 <C 1 and 
constitute small corrections to the standard ACDM expansion. 

In the quasi-static limit, the mass of the scalar field is given approximately by 



M 2 (a) 



1 



(2.26) 



3/rr(-Rgr) 

and hence we could in principle write the additional terms in the Friedmann equation solely 
in terms of M(a). However, this procedure is model dependent, as it requires writing / 
and /r in terms of /rr. 

With respect to the perturbations, neglecting all terms of order H 2 /M 2 and taking 
k(j)/a, kip /a 3> p, ip, the equations reduce to 



i 



2k 1 



3 + 2K 2 ) 



, 2l , r ,(3 + 2K 2 \ 2 

k 4> = -4ttG I — 3j ^ 2 I a p m 5 m , 



0, 



(2.27) 
(2.28) 
(2.29) 



where K = k/aM(a). In the quasi-static limit, terms of order H 2 /M 2 are small, so we have 
taken w = —1. We are considering modes that satisfy k 3> aH in equations (2.27 — 2.29); 
superhorizon modes will evolve according to General Relativity as they remain below the 
scalaron mass scale at all times (fc super h <C aM(a) for all a.) 

To summarize, we have argued that the quasi-static approximation, valid whenever 
B 2 I'M 2 <C 1, can be used to write the Friedmann equation in its standard General Rel- 
ativistic form with corrections of order 0(H 2 /M 2 ). In addition, all modifications to the 
perturbation equations can be expressed solely in terms of the scalar field mass M(a). In 
a forthcoming paper, the authors will consider how one might constrain the mass of the 
scalaron with cosmological data, using the above system of equations. 

2.4 The growth parameter in the quasi static approximation 

Now that we have a dynamical equation for 5 m , we can calculate 7. Performing the same 
steps as in the General Relativistic case, we find that 7 is described by the following 
equation 



Y(a,fc) = 



1 



ffii log[O m 

j --ill 



2 a 



+ 



2 2 

k 2 



2alog[n m ] \k 2 + a 2 M 2 



a 



a 



7 



log[^r 



(Am " 1) 



(2.30) 



This expression reduces to the General Relativistic equation (2.4) (with w = — 1) for 
M(a) —7- 00, as expected. The effect of modifying gravity is to introduce the last term on 
the right hand side of (2.30). 
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To obtain a solution valid for z > 1 we expand Q m ~ 1 — Q\oa 3 /Q m o + 0(a 6 ), in which 
case equation (2.30) reduces to 

7 ' ( °' fc) = "2a 7(a ' fc) + a " 2^ + aW)' ^ 

which has solution 

* i »4 + »'*-" ,! -^ , " ! (^5' < 2 - 32 ' 

where wo = Q,\o/£l m o and 70 is an arbitrary function of k, which is set to zero in order to 
recover the General Relativistic limit as a — > 0. We note that in deriving (2.31), we did 
not expand the term k 2 / (k 2 + a 2 M 2 (a)) in a <C 1, as it is a non-linear function of the scale 
factor 4 . 

For generic f(R) functions, the integral in (2.32) does not admit an analytic solution, 
however from (2.31) we can broadly describe the evolution of 7 for this class of models. 
At early times z > 10 the mass of the Scalaron will satisfy aM(a) » k for all modes of 
interest, and deviations from General Relativity will be suppressed. However, whenever 
any given k mode crosses the scalaron horizon k = a c M(a c ), the last term on the right 
hand side of (2.31) will dominate and 7 will decrease. There will be a zero in 7' due to the 
fact that the first and second terms grow relative to the third, and 7 will increase to the 
present, with asymptotic behaviour 7—^2/3 for a»l (this can be derived from the full 
equation (2.30).) 

Equation (2.31) and approximate solution (2.32) have been constructed by expanding 
the full equation (2.30) for a <C 1. Whilst this power series is strongly convergent and hence 
valid up to a ~ 1/2, we cannot obtain a solution in the region of particular interest (that 
is, z < 1) using this approach. However, as we shall see in section 3 our approximation is 
sufficient to accurately calculate the matter power spectrum at all redshifts. 

2.5 Application to specific f(R) models 

Now that we have an approximate, model independent form for 7 we can apply it to some 
existing f(R) functions in the literature. Specifically, we consider the functional forms 
(1.5,1.6). Our approach will be to solve the perturbation equation (2.29) numerically 
for <5 m , and use this in the definition (2.1) to obtain 7. We then compare this to our 
approximate 7, obtained by numerically integrating (2.32). 

Taking an initial redshift z\ = 100 and using model parameters n = 3/2, A = 5 x 10~ 4 , 
b = 2, g = 0.45, the results are exhibited in Figs. 1,2. We have neglected radiation, 
anticipating the fact that any deviations from General Relativity will occur for z < 10. We 
observe a close agreement between analytic and numerical solutions for z > 1, and a loss 
of accuracy in our approximation for z < 1, as expected. 

The dot-dashed curve in Figs. 1,2 is an approximate analytic solution for 7 in the a ~ 1 
region, which is discussed in the Appendix. For the power law model, we observe a close 



4 See [47] for earlier work on integral representations of 7 in modified gravity models. 
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agreement between the numerical and analytic solutions at low redshift, however for the 
exponential model our approximation for z < 1 is not particularly accurate. This is due to 
the fact that k modes of interest to the linear matter power spectrum only deviate from 
General Relativity at low redshift z < 1 for this model, whereas the key assumption made 
in obtaining this curve is that modifications to GR occur at z > 1. We direct the reader 
to the Appendix for further details. 

In Figs. 3, 4 we show j(a,k) for three different modes; k = 0.15, 0.1, 0.05/iMpc" 1 . For 
the model parameters chosen, 7 exhibits only a very weak k dependence at late times, 
taking the value 7 ~ 0.43. This value has been obtained previously in the literature for 
a number of different models (see for example [48]), and we derive it explicitly in the 
Appendix. 



0.6 




1+z 



Figure 1: The growth parameter 7 for the model (1.6), with model parameters n = 3, A = 5 x 10~ 4 
(solid), where we have taken k = 0.14/iMpc _1 . We have also exhibited approximate expressions for 
7 derived in section 2 in the a < 1 limit (dashed) and a ~ 1 (dot-dashed) limits. We observe a 
close agreement at all redshifts. 



3. Power spectrum reconstruction from 7 



Whilst our approximate form for 7, equation (2.32), is only accurate to ~ 10% at low 
redshift, we now use it to obtain the power spectrum P(a = l,k), and compare our re- 
sult to the power spectrum obtained by numerically evolving the full perturbation equa- 
tions (2.27 — 2.29). Taking the power law f(R) function (1.6) and the model parameters 
n = 1, A = 10~ 3 , our results are exhibited in Figs. 5, 6. In fig. 5 we exhibit the differ- 
ence #approx( a = l;^)/£m( a = 1j — 1> where 5 m is the density perturbation obtained by 
numerically evolving (2.29) and <5 a p pr ox is obtained by integrating the expression 



-'approx 



exp 



(3.1) 
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0.60 r 




_L , , 1 , 1 , , , , 1 , , , , 1 , , , , 1 1+z 

1 2 3 4 5 

Figure 2: The growth parameter 7 for the model (1.5), with model parameters b = 2, g = 0.45 
(solid), where we have taken k = 0.14/iMpc -1 . We have also exhibited approximate expressions for 
7 derived in section 2 in the a <C 1 limit (dashed) and a ~ 1 (dot-dashed) limits. 
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Figure 3: The growth parameter 7(0, k) for the model (1.6), taking model parameters as in figs.l 
for k = 0.15/iMpc- 1 (solid), k = 0.1/iMpc" 1 (dashed), k = 0.05/iMpc- 1 . (dot-dashed). 

using our approximate solution (2.32) for 7. We observe a close agreement between the 
two approaches, which suggests that our approximation (2.32) can be used to obtain the 
power spectrum for these models at any redshift, despite the apparent loss of accuracy at 
z ~ 0. We have confirmed this statement for a wide range of parameter choices. 

In fig. 6 we exhibit the power spectrum for the f(R) model (1.6) at z = 0. This has been 
obtained by evolving (2.29) in the range z = (20,0) to obtain 5 m (a,k), and normalizing 
the power spectrum to its General Relativistic value at z = 20. The resulting f(R) and 
General Relativistic power spectra are exhibited as solid and dotted lines respectively. We 
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Figure 4: The growth parameter 7(0, k) for the model (1.5), taking model parameters as in figs.l 
for k = 0.15/iMpc- 1 (solid), k = 0.1/iMpc" 1 (dashed), k = 0.05/iMpc^ 1 . (dot-dashed). 

observe an increased power on small scales for the modified gravity model. This is due to 
the fact that for the model parameters chosen, high k modes will satisfy k > aM(a) at 
late times. Whenever this occurs, there will be an order ~ O(10%) increase in G e s in the 
S m dynamical equation, leading to enhanced power. Low k modes, which do not satisfy 
k > aM{a) at any redshift, will not observe any increase in G e g and will evolve according 
to General Relativity. 

In fig. 6 we have exhibited the f(R) power spectrum using both the full equation (2.29) 
to obtain <5 m and the approximate solution <5 a ppro X ; the approximate solution is the dashed 
curve which closely overlaps the full solution in fig. 6. We observe an excellent agreement 
between the two methods. 

3.1 Comment on the 7 parameterization for f(R) gravity 

So far in this paper we have derived the expression (2.32) for the growth parameter 7, 
involving a model dependent integral over the mass of the scalaron. The M(a) dependence 
significantly alters the behaviour of 7 from its General Relativistic form; 7 will now be a 
non-analytic function of a, k, O to q, Hq and will depend on the explicit functional form of 
f(R). To obtain 5 m using 7, one must numerically integrate the expression (2.7) relating 7 
and 5 m , where 7 itself involves a non-analytic, model dependent integral over M(a). Given 
these facts, it is not obvious why one should persevere with 7 as a parameterization for 
growth in f(R) gravity. 

In addition, one should be careful not to construct an approximate functional form for 
7, valid over a particular redshift range. For example, in the previous section we found 
that 7 will exhibit only a very weak k dependence at low redshift. One might therefore 
be tempted to construct a simple functional form for 7 = 7(a) for a ~ 1 and use this in 
(2.7) to obtain 8 m . However, this approach will not uniquely specify the power spectrum, 
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Figure 5: The quantity (<5^ — <5approx)/^m> which is the fractional difference in 5^ using the two 
different approaches explained in the text. The fractional error is small (less than <~ 2%), suggesting 
that our expression (2.32) for 7 is suitable to reconstruct the power spectrum. 

as the k dependence of 5 m will remain undetermined. This is due to the fact that 7 is only 
denned via, the ratio S^fS-m (see (2.1)), and. so writing <5m cts $m = 5o(k)5i(a, k) (where 
S\(a, k) is only weakly k dependent at z ~ 0), it is clear that 7 will give us no information 
on 5o(k). In this respect, one can think of f(R) models as modifying the power spectrum 
at late times in two ways; changing the time evolution of the amplitude and introducing 
an additional component to the transfer function. Phenomenological constructions of 7(a) 
will yield information on the time evolution, but not the modified transfer function. Care 
should be taken for any model in which <5 m exhibits both scale and time dependence. 

4. Parameterizing the scalaron mass M(a) 

We have argued above that directly parameterizing 7(0, k) is problematic in f(R) gravity. 
Given that we can write the modified Friedmann and perturbation equations in terms of 
the scalaron mass M(a), it seems more appropriate to parameterize this quantity. Once 
we have a functional form for M(a), one can then use it in (2.32) to obtain an expression 
for 7, if one wishes to persevere with this parameterization. In doing so, both the a and k 
dependence of 7 and hence the power spectrum will be accounted for. 

The advantage of using M(a) is that it is not scale dependent, and the essential features 
of viable f(R) models can be captured with a very simple functional form. All we require 
is a function M(a) which monotonically grows to the past, such that i?GR/^ 2 ( a ) ~~ * 
as a —7- (one should introduce a cutoff to ensure that M{a) does not diverge as a — > 0, 
however this is unnecessary as we only consider the late universe z < 1000.) In fig. 7 we 
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Figure 6: Power spectrum for the model (1.6), with model parameters n — 1, A — 1 x 10~ 3 (solid 
line) and the General Relativistic power spectrum (dotted line). We observe an increased power 
at small scales, as expected. The dashed line, practically indistinguishable from the solid line, is 
the power spectrum obtained using our approximate scheme (that is, numerically integrating (2.32) 
and using the resulting 7 to obtain <5 ap prox from (3.1). 




Figure 7: The scalar field mass M(a)/H n for the exponential (dashed) and power law (solid) 
models (1.5, 1.6), with model parameters n — 3, A = 2 x 10~ 3 , g = 0.48, b = 3. Both monotonically 
increase to the past. 



show the simple behaviour of the mass M(cl)/Hq for the power law and exponential models, 
for parameter choices n = 3, A = 2 x 10 -3 , g = 0.48, 6 = 3 and O m o = 0.3. 

A simple functional form for M(a), which is representative of 'viable' f(R) models, is 
given by 
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M(a) = ^, (4.1) 

where jjl > 1 and a > 4 are free parameters. This mass corresponds to an f(R) model that 
is is free from instabilities for < z < 10 3 , and reduces to General Relativity at early times 
(subject to a > 4; this parameter dictates how quickly the mass diverges to the past.) 
It therefore constitutes a viable f(R) model in the regime of interest. We note that this 
parameterization has been considered previously; see [61]. 

Using this mass in (2.32), we can write the growth parameter 7 in terms of a hyper- 
geometric function 2-F1 [a, b, c; x], 



_ 6 k l n m0 

7 11 M 2 (a)a 5 O\ (l+4cT) 2 1 



4a + 1 8a -3 k 2 

J- 1 



(4.2) 



4fj-4'4f7-4' a 2 M 2 (a)_ ' 
Similarly, we can write the modified Friedmann equation (2.23) solely in terms of M(a), 



where the function a(a) is given by 

, \ 3Sl m o 
y ' 4a 6 



2a - 5 _ 4cr - 13 



2a -3 ' " m0 a -3 



(4.4) 



By using (4.2) in (3.1), numerically integrating this expression over low redshifts z < 20 and 
normalising the power spectrum to General Relativity at z = 20, one can obtain the linear 
matter power spectrum for this parameterization. Hence (4.2) and (4.3) are sufficient to 
fully parameterize the growth and expansion histories for the class of 'viable' f{R) models, 
subject to the quasi static approximation M{a) > H(a). 

To highlight the effect of M(a) on various observational probes, we exhibit the lumi- 
nosity distance, the CMB angular power spectrum and the matter power spectrum for the 
model (4.1) in figs. 8-10, taking model parameters ji = 6, a = 4. Specifically, in fig.8 we 
show the difference Am = 5(log[dL] — log[<i[ GR ' ) ]) between the luminosity distance for the 
modified gravity model (4.1) and its ACDM value d|° R \ where o?l is defined in the usual 
fashion 

*w-< 1+ "jC^- (4 ' 5) 

We observe practically no modified gravity signal in either the luminosity distance or the 
CMB angular power spectrum, indicating that the most stringent constraints on these 
models will be be derived from the matter power spectrum, as shown in Fig. 10. 



5. Discussion 

In this work we have considered the evolution of density perturbations for the class of 
so-called 'viable' f(R) modified gravity models. By using the quasi static approximation, 
we have written the Friedmann equation as an expansion around General Relativity, with 
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Figure 8: The difference between the luminosity distance for the modified gravity model character- 
ized by the mass (4.1) and the General Relativistic luminosity distance dj ^. There is practically 
no modified gravity signal in dh, indicating that the expansion history will be an ineffective probe 
of these models. 
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Figure 9: The CMB angular power spectrum for the model (4.1) (solid line) and the General 
Relativistic power spectrum (dashed line), taking Q mn — 0.29, Oao = 0.71, n s = 1, h = 0.7. We 
observe no significant modified gravity signature, even at large angular scales. 

corrections of order Hq^/M 2 . The perturbation equations can similarly be written in 
terms of the mass of the scalar field, and hence we can parameterize both the expansion 
and growth histories solely in terms of M(a). In a future publication, we will use the 
modified Friedmann and perturbation equations to derive constraints on the mass of the 
scalar field in f(R) models with cosmological probes. 

We have also constructed an approximate functional form for the growth parameter 
for this class of models, again in terms of the mass M{a). Whilst our approximate solution 
is only accurate to ~ 10% for < z < 1, it can still be used to accurately reconstruct 
the power spectrum at all redshifts. We have shown that by specifying a suitable mass 
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Figure 10: The matter power spectrum for the model (4.1) (solid line) and the General Relativistic 
power spectrum (dashed line), Sl m o = 0.29, f^Ao = 0.71, n s = 1, h = 0.7. We observe a significant 
increase in power at small scales, as expected. 

function M(o), equations (2.23) and (2.32) are sufficient to constrain this class of models 
with cosmological data. 

We would like to conclude by comparing our approach to existing parameterizations 
in the literature. The perturbation equations (2.27 — 2.29) belong to a parameterization 
considered in [59], [60], where arbitrary functions rj(a,k) and fJ,(a,k) were introduced into 
the perturbation equations 



Equations (5.1,5.2) have a broad range of applicability and represent an extremely gen- 
eral parameterization of modified gravity models, of which we have only considered a very 
specific subset. However, one advantage to our approach is that we have an action from 
which to derive the field equations, and hence both the background expansion and pertur- 
bation equations are consistent (it is not clear how the background expansion is modified 
if we simply choose rj(a,k) and fi(a,k).) We also note that when equations (5.1,5.2) are 
implemented, typically the k dependence of r\ and fi is neglected. 

Finally, we comment on a series of papers containing the most comprehensive cosmo- 
logical constraints on f(R) models to date [27], [28], [32]. In these works, the f(R) terms in 
the modified Friedmann equation (2.8) are written as the energy density of a fluid (typically 
a perfect fluid with constant equation of state), and this equation is treated as a second 
order differential equation for f(a). By demanding that the model reduces to General Rel- 
ativity as a —7- 0, one of the integration constants is fixed. There is only one remaining free 
parameter, taken to be Bq as defined in the introduction. Then the perturbation equations 
in all three domains of interest (superhorizon, subhorizon and non-linear) are written in 
terms of B(a), and Bq is constrained using cosmological data sets. This approach is more 





k 2 (f) = —4TrG/j,(a, fc)a 2 p m (5, 



(5.2) 
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sophisticated than ours, in the sense that the non-linear regime has also been considered 
(we will analyse the non-linear regime in a forthcoming publication, see also [62].) 

The most significant difference between this approach and ours is that we have not 
fixed the background expansion history; although it remains close to ACDM there will 
be corrections of order Hq R /M 2 . The mass of the scalar field M(a) therefore contains 
two parameters. One dictates the value of M(a) relative to the Hubble parameter at the 
present time (this is essentially Bq), and the other dictates how quickly M(a) grows to the 
past. By fixing the evolution of H(a), the time evolution of M(a) is fixed, and it can be 
shown that an f(R) model with equation of state w = — 1 corresponds to an f(R) model 
given by 



during the matter dominated epoch 1 < z < 10 3 , where A is the free model parameter 
(which is related to Bo), and p + = (—7 + \/73)/4 ~ 0.4. This is an f(R) model with 
particularly slow asymptotic behaviour to the past; it is of the form (1.6) with n ~ 0.06. 
It would be instructive to see what constraints on Bo can be obtained for models where 
M(a) possesses a steeper redshift dependence. This will be the subject of a forthcoming 
publication. 

Acknowledgments: The authors would like to thank Eric hinder and Shaun Thomas 
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6. Appendix: low redshift behaviour of 7 



To obtain an approximate expression for 7 in the region < z < 1, we use the simplifying 
assumption that 9 = {3+AK 2 ) / {3+3K 2 ) in (2.29) acts heaviside function, interpolating 
between 9(a,k) ~ 1 for k <C aM(a) and 6* ~ 4/3 for » aM(a). We therefore solve the 
equations 



a 2 6'^ ll (a, k) + 3a 
o 2 ^ 2 (a, k) + 3a 



1 - 
1 - 



^ nl (a, fc) = -ns ml , 

S' m2 (a, k) = 2JW m2 , 



(6.1) 
(6.2) 



and match e) m i,2 and S' ml 2 at a = a c (/c). If we assume that the relevant modes cross the 
scalaron horizon for z > 1, where the small a expansion can be used to match the two 
solutions at a = a c , we obtain the following solution for z < 1 



-So 



7 + V33 \ /fi 



'mO 



dcO 
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a i- m +(l-O m )(- 1+ ^/ 12 x 
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(6.3) 



12 



12 



, fi m [a] 



where the constants Ci and C2 are given by 
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and m + is given by 



m_|_ 



-1 + V33 



(6.4) 
(6.5) 
(6.6) 
(6.7) 
(6.8) 

(6.9) 



(6.10) 



We stress that (6.3) is only valid at late times, and only for k modes that have crossed 
the 'scalaron horizon' at early times z > 1. The key approximation that has been used in 
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deriving our solution; that 9(a, k) can be approximated by a discontinuous step function, 
is not applicable for modes that cross the horizon for z < 1. 

The solution (6.3) can be written as <5 m = 5oB(a)a c ~ m+ , where B(a) is a model 
independent function of a. All k and f(R) model dependence is incorporated in the al m+ 
term, which will introduce a tilt to the matter power spectrum, relative to the standard 
ACDM case. This tilt can be derived explicitly, by inverting the expression a c M{a c ) = k 
for any given model, and substituting the resulting a c (k) into our expression for <5 m . For the 
power law model (1.6) we find S m = AoB(a)k^"^~ 5 ^^ 12n+8 \ where Aq is an unimportant 
constant. This power law k dependence has been obtained in [50, 51, 48]; for the exponential 
model (1.5), <5 m acquires a logarithmic dependence. 

Our calculation is in agreement with existing work in the literature. Since 5 m is a 
separable function of k and a for z ~ 0, it follows that 7 is independent of k for z ~ 0, at 
the level of approximation to which we are working. Using (6.3) in the definition of 7 and 
taking Q m Q = 0.27, we obtain "f\ a =i — 0.42. Both of these results are in agreement with 
[48], where it was noted that at z = 0, 7 has a small dispersion in k and takes the value 
7 ~ 0.4 - 0.43. 
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